###################################################
#### "Grocery for Votes" ##########################
#### Replication File #############################
#### Randomization Inference (Figure 4) ###########
###################################################

rm(list=ls(all=TRUE))
library(foreign)
library(ggplot2)
library(stargazer)
library(arm)
library(rms)
library(msm)
library(systemfit)
library(xtable)
library(maptools)
library(ggplot2)
library(rgeos)
library(Imap)
library(spatstat)
library(multiwayvcov)
library(numDeriv)
library(interplot)


data<-read.dta("datawork.dta")
data<-data[data$distance<=20,]

#Finding centroids for each Polygon
centroids<-data[,c("centroidLONG","centroidLAT")]
coordinates(centroids) <- ~ centroidLONG + centroidLAT
centroids<-as.matrix(data.frame(centroids))


#Loading Soriana Information
sorianas<-read.dta('/rawdata/Coordenadas4Tiendas.dta')
coordinates(sorianas) <- ~ longitude + latitude
sorianas<-as.matrix(data.frame(sorianas))


placeboEPNslope<-c()
placeboEPNse<-c()
placeboAMLOslope<-c()
placeboAMLOse<-c()



for(k in 1:1000){
rm(a, min, distancesstores, benchmarkPRI.1, benchmarkPRD.1)  
a<-sorianas[sample(nrow(sorianas),71),]
min<-c()
for(i in 1:nrow(centroids)){
  distancesstores<-  matrix(gdist(centroids[i,1],centroids[i,2], a[,2],a[,1], units = "km"))
  min[i]<-min(distancesstores, na.rm=TRUE)
}
  

data$placebodist<-min
data$placeboprox<-1/data$placebodist

data1<-data
       
  
    #EPN
benchmarkPRI.1 <- lm( EPNa ~ placeboprox*PRIstr*highTURNOUT+placeboprox*PRDstr*highTURNOUT+placeboprox*PANstr*highTURNOUT+     
                        PRI09a+PRD09a+PAN09a+PART09+
                        lnpop+P18+P65+density+area+INDIGENOUS+CATHOLIC+NONRELIGIOUS+
                        EDUCATION+POSTEDUC+ILLITERACY+PROM_HNV+
                        PEAprop+PEAfemale+NOINSURANCE+FEMALEJEFA+
                        PERROOM+DIRTFLOOR+SERVICES+NO_SERVICES+
                        CAR+CELULAR+INTERNET+
                        id_mun, data=data1, x=TRUE, y=TRUE)

covBMpri <- cluster.vcov(benchmarkPRI.1, data1$id_dist)

coeftest(benchmarkPRI.1, covBMpri)

k1<-c()
k2<-c()
k3<-c()
k4<-c()

for(i in 1:length(benchmarkPRI.1$coef)){
  if(names(benchmarkPRI.1$coef[i])=="placeboprox"){
    k1<-i
  }
  if(names(benchmarkPRI.1$coef[i])=="placeboprox:highTURNOUT"){
    k2<-i
  }
  if(names(benchmarkPRI.1$coef[i])=="placeboprox:PRDstr"){
    k3<-i
  }
  if(names(benchmarkPRI.1$coef[i])=="placeboprox:highTURNOUT:PRDstr"){
    k4<-i
  }
}



g_prd <- function(b){
  return( b[k1] + b[k2] + b[k3] + b[k4] )
}

grad_g_prd <-  jacobian(g_prd, benchmarkPRI.1$coef)



placeboEPNslope[k] <- benchmarkPRI.1$coefficients[c("placeboprox")] + 
  benchmarkPRI.1$coefficients[c("placeboprox:highTURNOUT")] + 
  benchmarkPRI.1$coefficients[c("placeboprox:PRDstr")]+  
  benchmarkPRI.1$coefficients[c("placeboprox:highTURNOUT:PRDstr")]

placeboEPNse[k] <- sqrt(grad_g_prd%*% covBMpri %*% t(grad_g_prd))


#AMLO
    benchmarkPRD.1 <- lm( AMLOa ~ placeboprox*PRIstr*highTURNOUT+placeboprox*PRDstr*highTURNOUT+placeboprox*PANstr*highTURNOUT+     
                             PRI09a+PRD09a+PAN09a+PART09+
                             lnpop+P18+P65+density+area+INDIGENOUS+CATHOLIC+NONRELIGIOUS+
                             EDUCATION+POSTEDUC+ILLITERACY+PROM_HNV+
                             PEAprop+PEAfemale+NOINSURANCE+FEMALEJEFA+
                             PERROOM+DIRTFLOOR+SERVICES+NO_SERVICES+
                             CAR+CELULAR+INTERNET+
                             id_mun, data=data1, x=TRUE, y=TRUE)
    
    covBMprd <- cluster.vcov(benchmarkPRD.1, data1$id_dist)
    
    coeftest(benchmarkPRD.1, covBMprd)
    
    k1<-c()
    k2<-c()
    k3<-c()
    k4<-c()
    
    for(i in 1:length(benchmarkPRD.1$coef)){
      if(names(benchmarkPRD.1$coef[i])=="placeboprox"){
        k1<-i
      }
      if(names(benchmarkPRD.1$coef[i])=="placeboprox:highTURNOUT"){
        k2<-i
      }
      if(names(benchmarkPRD.1$coef[i])=="placeboprox:PRDstr"){
        k3<-i
      }
      if(names(benchmarkPRD.1$coef[i])=="placeboprox:highTURNOUT:PRDstr"){
        k4<-i
      }
    }
    
    
    
    g_prd <- function(b){
      return( b[k1] + b[k2] + b[k3] + b[k4] )
    }
    
    grad_g_prd <-  jacobian(g_prd, benchmarkPRD.1$coef)
    
    
    placeboAMLOslope[k] <- benchmarkPRD.1$coefficients[c("placeboprox")] + benchmarkPRD.1$coefficients[c("placeboprox:highTURNOUT")] + benchmarkPRD.1$coefficients[c("placeboprox:PRDstr")]+  benchmarkPRD.1$coefficients[c("placeboprox:highTURNOUT:PRDstr")]
    
    placeboAMLOse[k] <- sqrt(grad_g_prd%*% covBMprd %*% t(grad_g_prd))
    
print(k)
  
  }
  

save(placeboEPNslope, file="EPNslopeRI.RData")
save(placeboAMLOslope, file="AMLOslopeRI.RData")
save(placeboEPNse, file="EPNseRI.RData")
save(placeboAMLOse, file="AMLOseRI.RData")


benchmarkPRI.1 <- lm( EPNa ~ proximity*PRIstr*highTURNOUT+proximity*PRDstr*highTURNOUT+proximity*PANstr*highTURNOUT+     
                        PRI09a+PRD09a+PAN09a+PART09+
                        lnpop+P18+P65+density+area+INDIGENOUS+CATHOLIC+NONRELIGIOUS+
                        EDUCATION+POSTEDUC+ILLITERACY+PROM_HNV+
                        PEAprop+PEAfemale+NOINSURANCE+FEMALEJEFA+
                        PERROOM+DIRTFLOOR+SERVICES+NO_SERVICES+
                        CAR+CELULAR+INTERNET+
                        id_mun, data=data, x=TRUE, y=TRUE)

covBMpri <- cluster.vcov(benchmarkPRI.1, data$id_dist)

coeftest(benchmarkPRI.1, covBMpri)

g_prd <- function(b){
  return( b[2] + b[140] + b[142] + b[147] )
}

grad_g_prd <-  jacobian(g_prd, benchmarkPRI.1$coef)

slopesEPN <- benchmarkPRI.1$coefficients[c("proximity")] + benchmarkPRI.1$coefficients[c("proximity:highTURNOUT")] + benchmarkPRI.1$coefficients[c("proximity:PRDstr")]+  benchmarkPRI.1$coefficients[c("proximity:highTURNOUT:PRDstr")]

seEPN<-sqrt(grad_g_prd%*% covBMpri %*% t(grad_g_prd))

benchmarkPRD.1 <- lm( AMLOa ~ proximity*PRIstr*highTURNOUT+proximity*PRDstr*highTURNOUT+proximity*PANstr*highTURNOUT+     
                        PRI09a+PRD09a+PAN09a+PART09+
                        lnpop+P18+P65+area+density+INDIGENOUS+CATHOLIC+NONRELIGIOUS+
                        EDUCATION+POSTEDUC+ILLITERACY+PROM_HNV+
                        PEAprop+PEAfemale+NOINSURANCE+FEMALEJEFA+
                        PERROOM+DIRTFLOOR+SERVICES+NO_SERVICES+
                        CAR+CELULAR+INTERNET+
                        id_mun, data=data, x=TRUE, y=TRUE)

covBMprd <- cluster.vcov(benchmarkPRD.1, data$id_dist)

coeftest(benchmarkPRD.1, covBMprd)

g_prd <- function(b){
  return( b[2] + b[140] + b[142] + b[147] )
}

grad_g_prd <-  jacobian(g_prd, benchmarkPRD.1$coef)



slopesAMLO <- benchmarkPRD.1$coefficients[c("proximity")] + benchmarkPRD.1$coefficients[c("proximity:highTURNOUT")] + benchmarkPRD.1$coefficients[c("proximity:PRDstr")]+  benchmarkPRD.1$coefficients[c("proximity:highTURNOUT:PRDstr")]

seAMLO<-sqrt(grad_g_prd%*% covBMprd %*% t(grad_g_prd))

placeboEPN<-placeboEPNslope
placeboAMLO<-placeboAMLOslope

EPN<-slopesEPN
AMLO<-slopesAMLO

placebo<-c(placeboEPNslope,placeboAMLOslope)
placebo<-data.frame(placebo)
placebo$reg<-NA
placebo$reg[1:1000]<-"Peña Nieto's Vote Share"
placebo$reg[1001:2000]<-"López Obrador's Vote Share"
placebo$val<-NA
placebo$val[1:1000]<-slopesEPN
placebo$val[1001:2000]<-slopesAMLO

p    <- 0.05


#Randomization Inference Results
placebo$pvalue<-NA
placebo$pvalue[1:1000]<- mean(placebo$placebo[1:1000] > EPN)  
placebo$pvalue[1001:2000]<-  mean(placebo$placebo[1001:2000] < AMLO)  


placebo$val1[1:1000]<-slopesEPN
placebo$val1[1001:2000]<-slopesAMLO

placebo$placebo1[1:1000]<-placeboEPNslope
placebo$placebo1[1001:2000]<-placeboAMLOslope


placebo$pvalue1<-NA
placebo$pvalue1[1:1000]<- mean(placebo$placebo1[1:1000] > slopesEPN)  
placebo$pvalue1[1001:2000]<-  mean(placebo$placebo1[1001:2000] < slopesAMLO)  



placebo$adjust<-NA
placebo$adjust[1:1000]<- (7.5) 
placebo$adjust[1001:2000]<-(-7.5)




p <- ggplot(placebo, aes(x=placebo))+geom_histogram(binwidth=.5, colour="black", fill="white")+
  theme_bw()+geom_vline(aes(xintercept=val), linetype=2, size=1, colour="red")+ylab("Density")+
  xlab(expression(paste(partialdiff, "y", "/" , partialdiff, "Proximity")))+ facet_grid(reg ~ .)+
  theme(axis.text=element_text(size=16),
  axis.title=element_text(size=16),
  strip.text = element_text(size = 16))


